function eigenvalue_experiment()
% EIGENVALUE_EXPERIMENT 特征值算法综合实验
% 
% 本实验涵盖特征值计算的各种算法和应用

fprintf('=== 特征值算法综合实验 ===\n\n');

%% 实验1：幂迭代法基础
fprintf('实验1：幂迭代法基础\n');
fprintf('目标：理解幂迭代法的基本原理和收敛性\n\n');

% 测试矩阵
A1 = [4, 1, 0; 1, 3, 1; 0, 1, 2];
fprintf('测试矩阵 A1 =\n');
disp(A1);

% 精确特征值
eigenvals_exact = sort(eig(A1), 'descend');
fprintf('精确特征值: [%.6f, %.6f, %.6f]\n', eigenvals_exact);

% 幂迭代实现
x0 = [1; 1; 1];  % 初始向量
max_iter = 20;
tol = 1e-10;

fprintf('\n幂迭代过程:\n');
fprintf('迭代    特征值近似    误差        收敛率\n');
fprintf('----    ----------    ----        ------\n');

x = x0 / norm(x0);
lambda_old = 0;
errors = [];

for k = 1:max_iter
    y = A1 * x;
    lambda = x' * y;  % Rayleigh商
    x = y / norm(y);  % 归一化
    
    error = abs(lambda - eigenvals_exact(1));
    errors = [errors; error];
    
    % 计算收敛率
    if k > 1
        conv_rate = errors(k) / errors(k-1);
    else
        conv_rate = NaN;
    end
    
    if k <= 10 || error < tol
        fprintf('%4d    %10.6f    %.2e    %8.4f\n', k, lambda, error, conv_rate);
    end
    
    if abs(lambda - lambda_old) < tol
        fprintf('收敛于迭代 %d\n', k);
        break;
    end
    
    lambda_old = lambda;
end

% 理论收敛率
theoretical_rate = abs(eigenvals_exact(2) / eigenvals_exact(1));
fprintf('理论收敛率: %.6f\n', theoretical_rate);
if k > 5
    actual_rate = mean(errors(6:end) ./ errors(5:end-1));
    fprintf('实际平均收敛率: %.6f\n', actual_rate);
end

%% 实验2：反幂迭代法
fprintf('\n实验2：反幂迭代法\n');
fprintf('目标：计算最小特征值和指定特征值\n\n');

% 计算最小特征值
fprintf('2.1 计算最小特征值:\n');
x_inv = x0 / norm(x0);

for k = 1:10
    y = A1 \ x_inv;  % 求解线性系统
    mu = x_inv' * y;
    x_inv = y / norm(y);
    
    lambda_inv = 1/mu;
    error_inv = abs(lambda_inv - eigenvals_exact(3));
    
    fprintf('迭代 %2d: λ = %8.6f, 误差 = %.2e\n', k, lambda_inv, error_inv);
    
    if error_inv < 1e-10
        break;
    end
end

% 带位移的反幂迭代
fprintf('\n2.2 带位移的反幂迭代 (目标: λ ≈ 3.2):\n');
target = 3.2;
sigma = target;
A_shift = A1 - sigma * eye(size(A1));

x_shift = x0 / norm(x0);
for k = 1:8
    y = A_shift \ x_shift;
    mu = x_shift' * y;
    x_shift = y / norm(y);
    
    lambda_shift = 1/mu + sigma;
    [~, closest_idx] = min(abs(eigenvals_exact - lambda_shift));
    error_shift = abs(lambda_shift - eigenvals_exact(closest_idx));
    
    fprintf('迭代 %2d: λ = %8.6f, 误差 = %.2e\n', k, lambda_shift, error_shift);
end

%% 实验3：Rayleigh商迭代
fprintf('\n实验3：Rayleigh商迭代\n');
fprintf('目标：实现三次收敛的特征值算法\n\n');

x_ray = x0 / norm(x0);
lambda_ray = x_ray' * A1 * x_ray;

fprintf('Rayleigh商迭代 (三次收敛):\n');
fprintf('迭代    特征值近似    误差        收敛率\n');
fprintf('----    ----------    ----        ------\n');

ray_errors = [];
for k = 1:6
    A_ray = A1 - lambda_ray * eye(size(A1));
    
    try
        y = A_ray \ x_ray;
        x_ray = y / norm(y);
        lambda_ray = x_ray' * A1 * x_ray;
        
        [~, closest_idx] = min(abs(eigenvals_exact - lambda_ray));
        error_ray = abs(lambda_ray - eigenvals_exact(closest_idx));
        ray_errors = [ray_errors; error_ray];
        
        if k > 1 && ray_errors(k-1) > 0
            conv_rate = ray_errors(k) / ray_errors(k-1)^3;  % 三次收敛
        else
            conv_rate = NaN;
        end
        
        fprintf('%4d    %10.6f    %.2e    %8.2e\n', k, lambda_ray, error_ray, conv_rate);
        
    catch
        fprintf('迭代 %d: 矩阵接近奇异\n', k);
        break;
    end
end

%% 实验4：QR迭代法基础
fprintf('\n实验4：QR迭代法基础\n');
fprintf('目标：同时计算所有特征值\n\n');

A4 = A1;
fprintf('基本QR迭代:\n');
fprintf('迭代    对角元素                    收敛指标\n');
fprintf('----    --------                    --------\n');

for k = 1:30
    [Q, R] = qr(A4);
    A4 = R * Q;
    
    convergence_measure = norm(tril(A4, -1), 'fro');
    
    if k <= 8 || mod(k, 5) == 0 || convergence_measure < 1e-10
        diag_elements = diag(A4);
        fprintf('%4d    [%6.4f, %6.4f, %6.4f]    %.2e\n', ...
                k, diag_elements, convergence_measure);
    end
    
    if convergence_measure < 1e-12
        fprintf('QR迭代收敛于第 %d 次迭代\n', k);
        break;
    end
end

final_eigenvals = sort(diag(A4), 'descend');
fprintf('QR迭代结果: [%.6f, %.6f, %.6f]\n', final_eigenvals);
fprintf('与精确值的误差: [%.2e, %.2e, %.2e]\n', ...
        abs(final_eigenvals - eigenvals_exact));

%% 实验5：带位移的QR迭代
fprintf('\n实验5：带位移的QR迭代\n');
fprintf('目标：加速QR迭代的收敛\n\n');

A5 = A1;
fprintf('带位移QR迭代:\n');

for k = 1:15
    n = size(A5, 1);
    shift = A5(n, n);  % 使用右下角元素作为位移
    
    [Q, R] = qr(A5 - shift * eye(n));
    A5 = R * Q + shift * eye(n);
    
    convergence_shift = norm(tril(A5, -1), 'fro');
    
    if k <= 8 || convergence_shift < 1e-10
        fprintf('迭代 %2d: 位移 = %6.4f, 收敛指标 = %.2e\n', ...
                k, shift, convergence_shift);
    end
    
    if convergence_shift < 1e-12
        fprintf('带位移QR收敛于第 %d 次迭代\n', k);
        break;
    end
end

%% 实验6：对称矩阵的Jacobi方法
fprintf('\n实验6：对称矩阵的Jacobi方法\n');
fprintf('目标：使用Jacobi旋转计算对称矩阵的特征值\n\n');

% 构造对称矩阵
A6 = [4, 1, 0; 1, 3, 1; 0, 1, 2];
fprintf('对称矩阵 A6 =\n');
disp(A6);

% Jacobi方法实现
V6 = eye(size(A6));  % 累积旋转矩阵
A_jacobi = A6;
max_jacobi_iter = 50;

fprintf('Jacobi方法迭代:\n');
fprintf('迭代    最大非对角元    收敛指标\n');
fprintf('----    ----------    --------\n');

for k = 1:max_jacobi_iter
    % 找到最大的非对角元素
    [max_val, max_idx] = max(abs(triu(A_jacobi, 1) + tril(A_jacobi, -1)), [], 'all', 'linear');
    [p, q] = ind2sub(size(A_jacobi), max_idx);
    
    if p > q
        temp = p; p = q; q = temp;  % 确保 p < q
    end
    
    if abs(A_jacobi(p, q)) < 1e-12
        fprintf('Jacobi方法收敛于第 %d 次迭代\n', k);
        break;
    end
    
    % 计算旋转角度
    if abs(A_jacobi(p, p) - A_jacobi(q, q)) < 1e-12
        theta = pi/4;
    else
        theta = 0.5 * atan(2 * A_jacobi(p, q) / (A_jacobi(p, p) - A_jacobi(q, q)));
    end
    
    c = cos(theta);
    s = sin(theta);
    
    % 构造Jacobi旋转矩阵
    J = eye(size(A_jacobi));
    J(p, p) = c; J(p, q) = -s;
    J(q, p) = s; J(q, q) = c;
    
    % 应用旋转
    A_jacobi = J' * A_jacobi * J;
    V6 = V6 * J;
    
    % 计算收敛指标
    off_diag_norm = norm(A_jacobi - diag(diag(A_jacobi)), 'fro');
    
    if k <= 10 || mod(k, 10) == 0 || off_diag_norm < 1e-10
        fprintf('%4d    %10.2e    %8.2e\n', k, max_val, off_diag_norm);
    end
end

jacobi_eigenvals = sort(diag(A_jacobi), 'descend');
fprintf('Jacobi方法结果: [%.6f, %.6f, %.6f]\n', jacobi_eigenvals);

% 验证特征向量
fprintf('特征向量验证 (||A*v - λ*v||):\n');
for i = 1:3
    v = V6(:, i);
    lambda = jacobi_eigenvals(i);
    residual = norm(A6 * v - lambda * v);
    fprintf('特征值 %d: %.2e\n', i, residual);
end

%% 实验7：大规模稀疏矩阵
fprintf('\n实验7：大规模稀疏矩阵特征值\n');
fprintf('目标：处理大规模稀疏矩阵的特征值问题\n\n');

% 构造稀疏三对角矩阵
n = 100;
e = ones(n, 1);
A7 = spdiags([e, 2*e, e], [-1, 0, 1], n, n);
fprintf('稀疏三对角矩阵: %d×%d\n', n, n);
fprintf('非零元素数: %d (密度: %.2f%%)\n', nnz(A7), 100*nnz(A7)/n^2);

% 使用MATLAB的eigs函数计算部分特征值
k_eigs = 6;  % 计算前6个特征值
fprintf('计算前 %d 个特征值...\n', k_eigs);

tic;
[V_sparse, D_sparse] = eigs(A7, k_eigs, 'largestabs');
time_eigs = toc;

sparse_eigenvals = diag(D_sparse);
fprintf('稀疏特征值计算时间: %.4f 秒\n', time_eigs);
fprintf('前%d个特征值: [%.4f, %.4f, %.4f, %.4f, %.4f, %.4f]\n', ...
        k_eigs, sparse_eigenvals);

% 与完整特征值分解比较（小规模）
if n <= 50
    tic;
    full_eigenvals = eig(full(A7));
    time_full = toc;
    
    fprintf('完整特征值分解时间: %.4f 秒\n', time_full);
    fprintf('时间比: %.2f\n', time_full / time_eigs);
end

%% 实验8：条件数和扰动分析
fprintf('\n实验8：条件数和扰动分析\n');
fprintf('目标：分析特征值对矩阵扰动的敏感性\n\n');

% 良态矩阵
A8_good = diag([1, 2, 3, 4]);
cond_good = cond(A8_good);
fprintf('良态对角矩阵条件数: %.2f\n', cond_good);

% 病态矩阵
epsilon = 1e-8;
A8_ill = [1, 1; epsilon, 1];
cond_ill = cond(A8_ill);
fprintf('病态矩阵条件数: %.2e\n', cond_ill);

% 扰动分析
perturbation_levels = [1e-10, 1e-8, 1e-6, 1e-4];
fprintf('\n扰动分析结果:\n');
fprintf('扰动水平    特征值变化    相对变化\n');
fprintf('--------    ----------    --------\n');

eigenvals_original = eig(A8_ill);
for pert_level = perturbation_levels
    E = pert_level * randn(size(A8_ill));
    A_perturbed = A8_ill + E;
    eigenvals_perturbed = eig(A_perturbed);
    
    % 计算特征值变化
    eigenval_change = norm(sort(eigenvals_perturbed) - sort(eigenvals_original));
    relative_change = eigenval_change / norm(eigenvals_original);
    
    fprintf('%8.0e    %10.2e    %8.2e\n', pert_level, eigenval_change, relative_change);
end

%% 实验9：应用实例 - 主成分分析
fprintf('\n实验9：应用实例 - 主成分分析\n');
fprintf('目标：使用特征值分解进行数据降维\n\n');

% 生成相关数据
n_samples = 100;
X_raw = randn(n_samples, 3);
% 创建相关性
X_raw(:, 2) = X_raw(:, 1) + 0.5*X_raw(:, 2);
X_raw(:, 3) = 0.3*X_raw(:, 1) + 0.2*X_raw(:, 2) + 0.5*X_raw(:, 3);

fprintf('数据维度: %d×%d\n', size(X_raw));

% 中心化数据
X_centered = X_raw - mean(X_raw);

% 计算协方差矩阵
C = (X_centered' * X_centered) / (n_samples - 1);
fprintf('协方差矩阵:\n');
disp(C);

% 特征值分解
[V_pca, D_pca] = eig(C);
[eigenvals_pca, idx] = sort(diag(D_pca), 'descend');
V_pca = V_pca(:, idx);

fprintf('主成分特征值: [%.4f, %.4f, %.4f]\n', eigenvals_pca);
variance_explained = eigenvals_pca / sum(eigenvals_pca) * 100;
fprintf('方差解释比例: [%.1f%%, %.1f%%, %.1f%%]\n', variance_explained);

% 降维到2D
X_pca = X_centered * V_pca(:, 1:2);
fprintf('降维后数据维度: %d×%d\n', size(X_pca));

% 重构误差
X_reconstructed = X_pca * V_pca(:, 1:2)' + mean(X_raw);
reconstruction_error = norm(X_raw - X_reconstructed, 'fro') / norm(X_raw, 'fro');
fprintf('2D重构相对误差: %.4f\n', reconstruction_error);

%% 实验10：性能比较
fprintf('\n实验10：性能比较\n');
fprintf('目标：比较不同特征值算法的性能\n\n');

test_sizes = [10, 20, 50];
fprintf('矩阵规模    幂迭代    QR迭代    Jacobi    MATLAB eig\n');
fprintf('--------    ------    ------    ------    ----------\n');

for i = 1:length(test_sizes)
    n_test = test_sizes(i);
    A_test = randn(n_test, n_test);
    A_test = (A_test + A_test') / 2;  % 对称化
    
    % 幂迭代（只计算主特征值）
    tic;
    x_power = randn(n_test, 1);
    x_power = x_power / norm(x_power);
    for k = 1:100
        y = A_test * x_power;
        lambda_power = x_power' * y;
        x_power = y / norm(y);
    end
    time_power = toc;
    
    % QR迭代
    tic;
    A_qr_test = A_test;
    for k = 1:200
        [Q_test, R_test] = qr(A_qr_test);
        A_qr_test = R_test * Q_test;
        if norm(tril(A_qr_test, -1), 'fro') < 1e-8
            break;
        end
    end
    time_qr = toc;
    
    % Jacobi方法（对称矩阵）
    tic;
    V_jac_test = eye(n_test);
    A_jac_test = A_test;
    for k = 1:500
        [max_val, max_idx] = max(abs(triu(A_jac_test, 1) + tril(A_jac_test, -1)), [], 'all', 'linear');
        if max_val < 1e-8
            break;
        end
        [p, q] = ind2sub(size(A_jac_test), max_idx);
        if p > q
            temp = p; p = q; q = temp;
        end
        
        if abs(A_jac_test(p, p) - A_jac_test(q, q)) < 1e-12
            theta = pi/4;
        else
            theta = 0.5 * atan(2 * A_jac_test(p, q) / (A_jac_test(p, p) - A_jac_test(q, q)));
        end
        
        c = cos(theta); s = sin(theta);
        J_test = eye(n_test);
        J_test(p, p) = c; J_test(p, q) = -s;
        J_test(q, p) = s; J_test(q, q) = c;
        
        A_jac_test = J_test' * A_jac_test * J_test;
        V_jac_test = V_jac_test * J_test;
    end
    time_jacobi = toc;
    
    % MATLAB内置eig
    tic;
    eigenvals_matlab = eig(A_test);
    time_matlab = toc;
    
    fprintf('%6d      %6.4f    %6.4f    %6.4f    %10.4f\n', ...
            n_test, time_power, time_qr, time_jacobi, time_matlab);
end

%% 实验总结
fprintf('\n=== 实验总结 ===\n');
fprintf('1. 幂迭代法简单但只能计算主特征值\n');
fprintf('2. 反幂迭代和Rayleigh商迭代可以计算指定特征值\n');
fprintf('3. QR迭代法是计算所有特征值的通用方法\n');
fprintf('4. Jacobi方法特别适合对称矩阵\n');
fprintf('5. 对于大规模稀疏矩阵，使用专门的算法更高效\n');
fprintf('6. 特征值的条件数影响计算精度\n');
fprintf('7. 特征值分解在数据分析中有重要应用\n\n');

end